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Abstract 

Proton  exchange  membrane  fuel  cells  (PEMFCs)  have  attracted  much  attention  in  these  years.  In  PEMFCs,  liquid/gas  two-phase  flow  is  a 
common  phenomenon,  which  has  great  influence  on  fuel  cell  performance.  However,  the  liquid  water  transport  process  has  not  been  satisfactorily 
modeled  yet.  In  this  work,  a  two  dimensional  partial  flooding  model  was  developed,  in  which  the  pore  size  distribution  of  the  gas  diffusion  layer 
(GDL)  is  taken  into  consideration  in  the  explanation  of  fuel  cell  flooding  for  the  first  time.  Liquid  water  produced  is  considered  to  flood  a  fraction 
of  the  GDL  hydrophobic  pores  with  diameter  greater  than  the  capillary  condensation  threshold  diameter,  and  the  unflooded  pores  will  serve  as 
passageway  for  gas  transportation  and  the  corresponding  catalyst  area  is  available  for  electrochemical  reaction.  Use  this  model,  it  is  simple  to 
explain  membrane  dehydration  and  electrode  flooding.  Different  operation  conditions  have  been  studied  with  the  model  and  the  model  polarization 
curves  show  reasonable  accordance  with  the  experimental  results. 

©  2005  Elsevier  B.V.  All  rights  reserved. 
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1.  Introduction 

Proton  exchange  membrane  fuel  cell  (PEMFC),  named  by 
utilizing  a  proton  conducting  membrane  as  the  electrolyte  in  the 
electrode  sandwich,  is  one  of  the  promising  technologies  for 
future  electric  vehicles,  distributed  power  stations  and  portable 
power  sources  [  1  ] .  To  date,  most  widely  used  electrolyte  in  PEM¬ 
FCs  is  perfluorinated  membrane,  with  excellent  conductivity  and 
durability,  like  well  known  Nation®  membrane  from  Dupont. 
However,  as  water  is  the  proton  carrier  in  this  type  of  ionomer, 
which  limits  the  operation  temperature  of  fuel  cell  in  a  range 
where  liquid  water  is  stable,  liquid/gas  two-phase  phenomenon 
is  unavoidable  and  brings  difficulty  in  water  management.  On 
the  one  hand,  the  membrane  needs  to  be  well  hydrated  to  per¬ 
form  good  proton  conductivity,  which  requires  humidification 
of  reactants  or  taking  measures  to  retain  water  in  gas  diffusion 
layers  (GDLs)  or  membrane;  on  the  other  hand,  the  GDLs  need 
to  be  flooding-proof  to  ensure  free  accesses  for  gaseous  reactants 
transportation. 

In  the  past  years,  many  models  were  developed  to  study  the 
transport  and  reaction  phenomenon  in  fuel  cells  [2-8]  and  have 
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been  well  reviewed  by  Yao  et  al.  [9].  The  pioneering  models 
assumed  that  water  was  in  gaseous  state  [2 —4],  which  is  not  true 
for  high  current  density  or  heavy  humidification  operation.  In 
the  past  years,  liquid  water  transport  became  one  of  the  major 
concerns  in  fuel  cell  modeling  and  multi-phase  models  had  been 
developed  to  study  liquid/gas  transport  in  porous  GDL  and  flow 
channel  by  He  et  al.  [10],  Berning  and  Djilali  [11],  Kimble  and 
Vanderborgh  [12],  You  and  Liu  [13],  Wang  et  al.  [14].  Recently, 
there  were  some  reports  on  flooding  in  catalyst  layer  [15,16]. 
These  two-phase  flow  models  assumed  that  both  the  liquid  and 
gas  phases  were  continuous,  but  no  experimental  picture  con¬ 
firmed  this  assumption  [9].  Wang  [17]  classified  the  two-phase 
flow  and  transport  modeling  into  two  distinct  approaches;  con¬ 
tinuum  method  and  pore-scale  method,  most  of  the  existing 
models  belong  to  the  former  approach,  where  average  porous 
media  properties  such  as  porosity,  permeability  and  pore  diam¬ 
eter  were  applied  in  the  assumed  homogenous  diffusion  layer. 
However,  diffusivity  of  liquid  water  and  gas  varies  quite  much 
for  different  pore  properties  as  shown  in  [  1 8] .  Pore-scale  method 
for  liquid  transport  modeling  might  be  helpful  in  furthering  the 
understanding  of  transport  phenomena  in  fuel  cells. 

Till  now,  no  model  had  taken  the  pore  diameter  distribution 
into  consideration.  However,  GDL  pore  size  distribution  has 
greater  influence  on  mass  transport  than  total  porosity  in  PEMFC 
as  reported  by  Kong  et  al.  [19].  The  experimental  results  showed 
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Nomenclature 

aw  activity  of  water 

d  thickness  (cm) 

I)  diffusing  coefficient  (cm2  s-1) 

/  unflooded  pore  area  fraction 

F  Faraday  constant  (96485  C  mol- 1 ) 

i  average  current  density  (A  cm-2) 

j  local  current  density  (A  cm-2) 

j°  exchange  current  density  (A  cm-2) 

k  conversion  coefficient  (1/101325  atmPa-1) 

L  channel  length  (cm) 

Mm  equivalent  weight  of  dry  membrane  (g  mol-1) 

Mi  molar  flow  rate  of  species  i  (mol  s-1) 

«(]  electro-osmotic  drag  coefficient 

Nj  flux  of  species  i  in  y  direction  (mol  s-1  cm-2) 
p  gas  pressure  (atm) 

Psat  saturation  water  vapor  pressure  (atm) 

r  radius  of  pores  (cm) 

ro  average  pore  radius  (cm) 

R  gas  constant  (8.314  Jmol-1  K-1) 

Rm  membrane  resistance  (Q  cm) 

T  temperature  (K) 

V  voltage  (V) 

w  channel  width  (cm) 

Xj  molar  fraction  of  species 

Greek 

a  net  water  drag  coefficient  from  anode  to  cathode 

(H20/H+) 

aa  anode  transfer  coefficient 

ac  cathode  transfer  coefficient 

5  surface  tension  of  water  (N  cm- 1 ) 

e  porosity  of  GDL 

ri  overpotential  (V) 

X  water  content  in  membrane  (FI2O/SO3-) 

6  contact  angle  (°) 

pm  density  of  dry  membrane  (g  cm-3) 

crm  conductivity  of  membrane  (S  m-1) 

X  volume  fraction  of  hydrophilic  pores 

coefficient 

a>  coefficient  in  pore  distribution  function 

1 jr  pore  distribution  function 

Subscripts 
a  anode 

c  cathode 

h  hydrogen 

m  membrane 

n  nitrogen 

o  oxygen 

oc  open  circuit 

w  water 


that  macropores  (5-20  pim)  in  diffusion  layer  are  thought  to  pre¬ 
vent  water  flooding  of  electrode.  In  this  paper,  a  two  dimensional 
model  was  developed  to  account  for  the  partial  flooding  phe¬ 
nomenon  in  fuel  cells  with  GDL  pore  size  distribution  being 
considered  for  the  first  time. 


2.  Basic  idea  of  the  model 


Carbon  paper  and  carbon  cloth  are  often  applied  as  GDL  in 
PEMFC  after  hydrophobic  treatment  with  Teflon.  During  the 
treatment,  Teflon  is  not  possible  to  cover  the  carbon  fibers  uni¬ 
formly,  so  there  will  be  different  hydrophobic  pore  property 
in  GDL  because  of  the  hydrophobicity  difference  between  the 
two  materials.  In  this  model,  both  the  hydrophilic  pores  and  the 
hydrophobic  pores  are  taken  into  consideration.  For  hydrophilic 
pores,  capillary  condensation  of  liquid  water  will  occur  before 
saturation.  For  example,  water  will  condense  at  0.8  PsM(T)  in 
lOnm  diameter  pores  [20].  For  hydrophobic  pores,  on  the  con¬ 
trary,  condensation  will  occur  in  some  extent  of  oversatura¬ 
tion.  The  overpressure  of  water  condensation  in  pores  could  be 
expressed: 


28  cos  6 

A  p  =  k - 

r 


(1) 


Inside  the  fuel  cell,  even  though  the  gas  is  not  saturated, 
water  will  condense  in  hydrophilic  micropores  because  of  the 
negative  overpressure.  With  the  increase  of  water  vapor  pres¬ 
sure,  water  will  condense  in  smaller  pores  first  and  then  bigger 
ones.  When  the  gas  is  fully  saturated,  all  the  hydrophilic  pores 
are  flooded  with  liquid  water.  With  further  increase  of  water 
vapor  pressure,  water  will  condense  in  hydrophobic  pores  when 
water  vapor  pressure  exceeds  condensation  threshold  pressure 
in  hydrophobic  pores  with  some  diameter.  As  the  overpressure 
value  is  positive  and  greater  for  smaller  pores,  water  will  first 
condense  in  bigger  pores  and  then  smaller  pores  with  the  increase 
of  oversaturation.  So  corresponding  to  different  saturation  con¬ 
dition  along  the  flow  channel  in  fuel  cell,  there  will  be  different 
threshold  condensation  pore  diameter,  which  will  result  in  dif¬ 
ferent  local  unflooded  pore  fraction  and  different  local  active 
area  for  electrochemical  reaction.  If  we  assume  that  water  will 
not  condense  in  the  flow  channel  and  water  vapor  in  the  flow 
channel  is  in  equilibrium  with  the  liquid  water  in  GDL,  there 
will  be  water  form  and  transport  balance.  Along  the  channel, 
the  GDL  is  partially  flooded  to  different  extent;  however,  the 
fuel  cell  could  be  in  steady  operation  without  apparent  flood¬ 
ing  phenomenon  being  observed.  The  partial  flooding  process 
is  illustrated  in  Fig.  1. 


3.  Model  development  and  experimental 


3.1.  Model  assumptions 


In  this  model,  some  assumptions  are  applied: 


1.  The  GDL  is  composed  of  a  series  of  pores  with  differ¬ 
ent  diameters.  There  are  hydrophilic  pores  and  hydrophobic 
pores;  at  every  local  position,  the  hydrophilic  pores  and  the 
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- ►  x 

[]  Unfloodedhydrophilic  pores  Q  Flooded  hydrophilic  pores 
J  Unflooded  hydrophobicpores  ^  Flooded  hydrophobicpores 


Fig.  1.  Illustration  of  a  fuel  cell  with  partial  flooded  GDLs. 


hydrophobic  pores  have  the  same  property  of  pore  diameter 
distribution. 

2.  Water  cannot  condense  in  flow  channels. 

3.  Oxygen  and  hydrogen  will  not  dissolve  in  water. 

4.  Catalyst  layer  is  assumed  to  be  an  ultra-thin  layer. 

5.  Reactants  in  flow  channel  is  assumed  to  be  plug  flow. 


3.2.2.  Mass  transport  in  GDL 

According  to  Maxwell-Stefan  equation,  molar  fraction  gra¬ 
dient  of  hydrogen  in  the  anode  GDL  could  be  expressed: 

dXh  1 

.  =  -^(*hAw, a  -  Xw,dNh)  (4) 

dy  cDhw 

where  the  flux  of  hydrogen  and  water  are  related  to  the  local 
current  density  j  and  net  water  drag  coefficient  a: 


Nw,  = 


Ah  = 


2  F 


So  Eq.  (4)  changes  to  be: 

dXh  =  _L _ J_ 

dy  cDhw  IF 


[(2a  +  l)Xh  -  1] 


(5) 

(6) 


For  a  specific  position  along  the  flow  channel,  if  local  current 
density  is  fixed,  a  will  be  fixed  and  Eq.  (6)  could  be  integrated 
directly  to  be: 


Xh 


channel 

h 


2a  +  1 


exp 


j  y(2a  +  1) 
2  F  cDhw 


2a  +  1’  Xw’a  1 


(7) 


In  the  first  assumption,  the  same  pore  size  distribution  infor¬ 
mation  is  assumed  for  hydrophilic  and  hydrophobic  pores  at  any 
location  of  the  GDL.  The  pores  are  assumed  to  be  in  the  same 
length.  The  second  assumption  excludes  the  behavior  of  water 
draining  with  droplets,  which  is  difficult  to  express  in  mathe¬ 
matical  model.  We  assume  that  water  can  only  condense  in  the 
GDL  pores,  not  in  the  flow  channel,  and  that  water  is  drained 
out  of  the  fuel  cell  in  gaseous  state  in  oversaturation.  Hydro¬ 
gen  and  oxygen  are  assumed  to  be  indissolvable  in  water,  so 
the  gases  could  not  transport  to  the  catalyst  layer  through  the 
flooded  pores.  In  this  model,  as  the  main  concern  is  the  flooding 
of  GDL,  gas  species  concentration  in  the  flow  channel  is  not  so 
important,  so  plug  flow  is  assumed. 


3.2.  Model  equations 


3.2.1.  Species  inflow  channels 

As  the  flow  type  being  simplified  to  be  plug  flow  [3],  when  the 
fuel  cell  is  working  with  current  density  j,  the  gas  species  con¬ 
servation  equation  along  the  flow  channel  could  be  expressed: 

{fh  =  — 2;  §Wia  =  — 4 a;  anode  channel 

£o  =  —  1 ;  §w,c  =  2  +  4a;  £n  =  0;  (2) 

cathode  channel 


Here  anode  and  cathode  reactants  are  considered  to  be  in  co-flow 
pattern.  The  net  water  drag  coefficient  a  means  that  there  will 
be  a  mole  H2O  transport  from  anode  to  cathode  together  with 
1 .0  mole  H+  apparently. 

Molar  fraction  of  species  in  flow  channels  are: 


Xi  = 


Mi 

EM 


For  cathode  GDL,  molar  fraction  gradients  for  the  species 
are: 

dX/  _  V  (XiNj  ~  XjNi 

dy  V  c  Do 


/  No  \ 

/ 

Nn 

4/ 

\  N w,C  / 

V 

This  equation  could  be  rewritten  as: 

/  _  _  (4a  +  2)X0 


X 


d 

d^ 


w,c 


(8) 


{  N O  \ 
Nn 

\  N w,c  ) 

j 

(4a  +  2)Xn 

4  F 

(4a  +  2)X0  ( 

c^nw 

Xw,c  (  (4a  +  2)Xn 

V  C/)qw 

c^ow  c^nw  / 

(9) 


The  binary  diffusion  coefficients  in  the  former  equations  are 
related  with  temperature  and  can  be  found  in  [21]: 


c Dhw  =  3.68  x  10“5(T/307) 
cDon  =  8.40  x  10“6(77293) 
c  Dnw  =  1.01  x  10“5(T/293) 
cDow  =  1.01  x  10“5(77293) 


^.10 

mol  cm' 1  s“ 1 

1.06 

mol  cm-1  s_1 

(1.07 

mol  cm' 1  s'" 1 

!  1 .08 

mol  cm-1  s_1 

(10) 


The  binary  diffusion  coefficients  are  modified  to  effective 
diffusivity  with  Bruggeman  relation  as  used  in  many  models 
[2,3]: 

Df=DijS1-5 


(3) 


(ID 


1232 


Z  Liu  et  al.  /  Journal  of  Power  Sources  158  (2006)  1229-1239 


Saturated  water  vapor  pressure  at  temperature  T  could  be 
expressed  in  [3]: 

log  Psat  =  -19.0437  +  0.10847  -  2.1022  x  10“472 

+  1.4454  x  10“7r3  (12) 


Once  the  molar  fractions  of  the  species  are  known,  we  know 
the  saturation  extent  inside  GDL  and  the  threshold  flooding  pore 
radius  according  to  Eq.  (1)  for  anode  and  cathode  GDL. 


ra  =  k 


28  cos  0 


-..catalyst 
PX  w,a 


rc=k 


28  cos  6 


-..catalyst 
pX  w,C 


(13) 


As  discussed  in  Section  2,  hydrophilic  pores  with  radius 
smaller  than  this  threshold  radius  or  hydrophobic  pores  with 
radius  larger  than  the  threshold  radius  are  not  flooded  with  liquid 
water  and  available  for  gas  transportation,  and  the  corresponding 
catalyst  area  is  available  for  electrochemical  reactions,  so  if  the 
pore  radius  distribution  of  the  GDL  is  known,  the  unflooded  pore 
area  fraction,  also  active  catalyst  surface  area  fraction,  could  be 
calculated  with  the  threshold  radius: 

,philic  =  If  Mr)  dr  f phobic  =  So3  Mr)  dr 

/0°°  Mr)  dr  ’  '  a  /Q°°  Vhi(r)  dr  ’ 

rphilic  _  frc  Vfc(r)dr  phobic  _  f0L  ^c(r)dr 

~  If  Mr)  dr’  /c  Jf  MO  dr 

/a=/PhlllCX  +  /aph0blC(l-X), 

/c  =  /phmcX+/cph0biC(l-X)  (15) 


/Phlllc  and  f?hoblc  are  unflooded  pore  fraction  of  hydrophilic 
pores  and  hydrophobic  pores.  The  integrations  of  the  pore  size 
distribution  function  i//  results  in  pore  volume,  because  the  pores 
are  assumed  to  be  in  the  same  length,  the  volume  ratio  is  also 
unflooded  pore  area  fraction.  The  unflooded  pore  area  fraction 
is  the  sum  of  hydrophilic  pore  fraction  and  hydrophobic  pore 
fraction. 

Pore  size  distribution  of  GDL  was  shown  to  have  micropores 
and  macropores  [19,22]  with  back  layer  on  carbon  cloth  GDL, 
and  macropores  (5-20  |xm)  same  to  perform  flooding-proof 
function  [19].  Pore  size  distribution  in  Toray  carbon  paper  was 
measured  and  shown  to  have  dominant  macropores  in  5-100  |xm 
range  without  active  carbon  back  layer,  and  the  average  pore 
diameter  was  measured  to  be  23,  19  and  17  pirn  for  0,  16  and 
35%  PTLE  content  in  carbon  paper  [23].  In  our  previous  study, 
pore  size  distribution  of  Toray  carbon  paper  (TGP-H-060)  with 
active  carbon  back  layer  was  measured  and  results  showed  that 
most  of  the  pores  are  in  0-20  |xm  diameter  range,  and  the  average 
pore  diameters  are  10,  7  and  5  pi  in  for  24,  35  and  42%  PTLE  in 
carbon  paper,  respectively  [24,25].  The  GDL  pore  size  distribu¬ 
tion  is  close  to  normal  distribution.  In  this  model,  we  use  normal 
distribution  function  to  fit  the  experimental  pore  size  distribu¬ 
tion  with  24%  PTLE  in  carbon  paper,  average  pore  diameter  of 
which  is  10  pim.  The  pore  size  distribution  function  is: 


MO  =  Mr)  = 


1  e—  l/2a)2(r— r0)2 

y/2jta> 


Fig.  2.  Experimental  GDL  pore  size  distribution  from  reference  [24]  and  normal 
distribution  function. 


Lig.  2  shows  the  experimental  distribution  result  and  the  fit 
curve,  where  fq  =  10  pim  and  o>=  I .  With  Eqs.  ( 14)— (16),  for  a 
specified  pore  radius,  r,  the  unflooded  pore  area  fraction /a  and 
fc  could  be  calculated. 


3.2.3.  Water  transport  through  membrane 

There  are  three  mechanisms  for  water  transport  through  mem¬ 
brane  [26]:  electro-osmotic  drag  with  proton  transport,  back 
diffusion  by  water  concentration  gradient  and  convection  by 
pressure  difference.  In  this  model,  water  convection  is  not  taken 
into  account  because  balanced  pressures  are  applied.  So  water 
transport  equation  in  membrane  could  be  written  as  [3]: 


j  Pm  n  dX  j 

n  d - Dw —  =  a— 

F  Mm  dy  F 


(17) 


In  Eq.  (16)  the  electro-osmotic  drag  coefficient  and  diffusion 
coefficient  are  from  [3,27],  respectively: 


n  d 


2.5 

22 


X 


Dw  =  2.1  x  10  3X  exp 


(18) 

(19) 


At  the  interface  of  membrane  and  catalyst,  activity  of  water 
in  membrane  phase  is  equal  to  that  in  gaseous  phase,  which  is 
related  to  the  partial  pressure  of  water  in  gas: 

aw  =  Xw-L-  (20) 

"sat 

Water  content  in  membrane,  X,  which  is  defined  as  the  ratio 
between  the  number  of  water  molecules  to  the  number  of  charge 
sites  (S03_H+),  is  related  to  water  activity  in  membrane  accord¬ 
ing  to  the  relation  given  in  [3]: 


0.043  +  17.81aw  —  39.85n2  +  36.0 a3,  0  <  aw  <  1 
14.0  +  1.4(aw  —  1)  1  <  aw  <  3 


(16) 


(21) 
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3.2.4.  Electrochemistry 

Butler- Volmer  equations  are  utilized  to  calculate  current  den¬ 
sity  as  listed  in  Eq.  (22).  Only  the  unflooded  catalyst  surface  will 
contribute  to  the  electrochemical  reaction,  therefore,  the  func¬ 
tion  of  the  catalyst  active  area  fraction  appears  in  these  equations. 


7a  —  /a  7  a 


vcat 

•  =  f  -0 
jc  JcJc  ^.ref 


r  fa,F  \  (  (1  —  af)F  \1 

rl^H-expl — 

■  (acF  \  (  (1  -ac)F  \1 

exp  \  ~rt  r,c)  —  exp  \ - Uf  r,c)\ 


(22) 


Cathode  exchange  current  density  is  considered  to  vary  with 
temperature  with  a  relation  from  literature  [28],  where  the  acti¬ 
vation  energy  A E  has  a  value  of  27.7  kJ  mol-1 . 


7c, 2  =  7c,  t  exp 


A£x 

(  1 

Ml 

R 

KTi 

Tl). 

(23) 


Cathode  transfer  coefficient  also  varies  with  temperature 
according  to  the  following  relationship  [28]: 

ac  =  0.495  +  2.3  x  10“3(r  -  300)  (24) 


To  determine  the  ohmic  drop  of  the  membrane,  empirical 
relation  from  [3]  is  used  to  calculate  the  conductivity: 


<rm  =  (5.139A  —  3.26)  x  10  3  exp 


1268  x 


So  the  membrane  resistance  could  be  calculated  by  integra¬ 
tion  over  the  membrane  thickness: 
rdm  | 

Rm  =  /  — dy  (26) 

Jo  °m 

With  the  overpotentials  known,  the  cell  voltage  is  easy  to 
know  by  subtraction  the  overpotentials  from  the  open  circuit 
voltage: 

Vcell  =  ^oc  *7c  7^  m  (27) 

The  thermodynamic  open  circuit  voltage  Voc  is  a  function  of 
temperature  and  reactants  partial  pressure  with  Nernst  equation 
shown  below  [28] : 

Voc  =  1-229  -  8.456  x  10“4(T  -  298.15) 

+  4.31  x  10-5rin(/?H2/45)  (28) 


density  j  and  net  water  drag  coefficient  a  are  first  guessed,  and 
molar  fraction  of  the  species  in  anode  GDL  could  then  be  calcu¬ 
lated  algebraically  in  Eq.  (7)  and  numerically  with  fourth-order 
Runge-Kutta  method  for  cathode  GDL  in  Eq.  (9).  Threshold  pore 
radii  for  anode  and  cathode  ra  and  rc  could  then  be  calculated 
and  then  the  unflooded  pore  fractions /a  and/c.  Water  content  in 
membrane  X  is  calculated  by  integrating  Eq.  (17)  in  y  direction 
numerically  with  fourth-order  Runge-Kutta  method  and  gain  a 
new  net  water  drag  coefficient  a.  If  the  difference  between  new 
a  and  the  old  one  is  greater  than  a  specified  error  10-3,  an  iter¬ 
ation  loop  will  be  performed  with  secant  method.  With  a  proper 
a  found,  overpotentials  and  cell  voltage  are  calculated.  If  the 
cell  voltage  difference  between  specified  value  and  calculated 
value  is  higher  than  error,  another  iteration  loop  is  used  to  search 
for  a  correct  current  density  j  with  secant  method.  This  process 
is  performed  for  each  discretized  fraction  from  the  inlet  to  the 
outlet,  and  then  the  average  current  density  could  be  calculated 
by  integrating  the  current  density  distribution  along  the  flow 
channel. 

3.4.  Experimental 

Carbon  paper  (Toray  TGP-H-060)  was  immersed  in  30% 
PTFE  latex  (Dupont)  for  10  min  and  dried  in  air  and  then 
heat-treated  at  350  °C  for  1  h,  after  sprayed  active  carbon 
(VulcanXC-72,  Cabot  Corp.)  for  back  layer,  heat-treated 
again.  Pt/C  electrocatalyst  (Johnson-Matthey)  was  dispersed 
in  ethanol  solution  and  sprayed  to  the  surface  of  the  back  layer, 
with  Pt  loading  0.4  mg  cm-2.  Then  the  GDL  was  treated  in 
vacuum  dryer  at  140  °C  for  1  h.  The  treated  GDL  was  then  cut 
in  5  cm2  foursquare  pieces.  Pretreated  Nation  112  membrane 
(Dupont)  with  5vol.%  H2O2  (Beijing  Chemical)  solution  and 
0.5  M  H2SO4  solution  was  sandwiched  with  two  pieces  of 
GDL  and  hot-pressed  at  5  MPa  to  form  a  membrane  electrode 
assembly  (MEA).  Fuel  cell  performance  was  measured  with 
Arbin  fuel  cell  test  stand.  Hydrogen  and  air  are  humidified  with 
membrane  humidifier  from  Beijing  LN  Power  Sources  Co.  Ltd. 

The  former  mentioned  pore  size  distribution  measurement 
in  [24,25]  was  conducted  with  9500  Mercury  Porosimetry 
(Micromeritics  Instrument  Corp.,  USA).  The  sample  was  the 
carbon  paper  with  active  carbon  back  layer  prepared  in  the  pro¬ 
cedure  mentioned  above. 

4.  Results  and  discussion 


With  the  current  density  distribution  j  known,  average  current 
density  of  the  fuel  cell  at  a  specified  voltage  value  could  be 
calculated  by  integrating  j  along  the  whole  length  of  the  flow 
channel: 

1  fL 

l  =  T  7<»  d*  (29) 

E  Jo 

3.3.  Solution  method 

The  calculation  domain  is  discretized  in  x  and  y  direction 
and  the  equations  are  solved  numerically  at  a  specified  cell  volt¬ 
age.  At  a  discretized  fraction  along  the  flow  channel,  current 


4.1.  Base  case 

Parameters  used  in  this  model  are  listed  in  Table  1 .  Operation 
conditions  of  the  base  case  are  listed  in  Table  2.  The  experimental 
fuel  cell  is  a  5  cm2  single  cell  with  co-flow  serpentine  channels. 
The  width  of  the  channel  and  the  rib  is  1  mm,  and  channel  length 
is  25  cm.  In  the  model,  the  channel  width  is  set  to  be  2  mm, 
which  is  the  sum  of  channel  width  rib  width.  It  needs  be  pointed 
out  that,  in  the  along  the  channel  model,  the  information  in  the 
cross  channel  direction  could  not  be  presented;  there  is  some 
difference  between  the  2D  model  and  3D  model,  which  will  be 
further  discussed  in  the  following  sections.  As  the  active  area  is 
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Table  1 

Parameters  used  in  the  model 


Parameters 

Symbol 

Value 

Unit 

Channel  width 

w 

0.2 

cm 

Channel  length 

L 

25 

cm 

GDL  thickness 

daDL 

200 

(jim 

Membrane  thickness 

dm 

50 

(jim 

Membrane  equivalent  weight 

Mm 

1100 

g  mol- 1 

Dry  membrane  density 

P 

1.0 

gem-3 

Water  surface  tension 

S 

7.28  x  lO^6 

J  cm-2 

Anode  transfer  coefficient 

aa 

0.5 

- 

Cathode  transfer  coefficient 

ac 

0.5 

- 

Anode  exchange  current  density 

j° 

1.0 

A  cm-2 

Cathode  exchange  current  density  (343  K) 

J°c 

0.002 

A  cm-2 

Table  2 

Operation  conditions  for  the  base  case 

Parameters 

Value 

Unit 

Temperature  of  fuel  cell 

70 

°C 

Relative  humidity  of  reactants 

90 

% 

Pressure  of  reactants 

0.2 

MPa 

Flux  of  TL 

0.10 

SLmin-1 

Flux  of  air 

0.20 

SLmin-1 

Average  pore  diameter  of  GDL 

10 

[xm 

Porosity  of  GDL 

60 

% 

Hydrophilic  pore  volume  fraction 

30 

% 

Contact  angle  of  hydrophilic  pores 

75 

° 

Contact  angle  of  hydrophobic  pores 

125 

° 

the  same  in  the  model  and  in  experimental  cell,  the  same  base 
case  conditions  are  applied  for  both  the  experimental  and  model. 

Fig.  3  shows  the  current  density  distribution  along  the  flow 
channel  in  base  case  at  different  discharge  voltages:  0.7, 0.6,  0.5 
and  0.4  V.  It  is  shown  that  the  current  density  is  relatively  uni¬ 
form  distributed  at  high  cell  voltage,  e.g.  0.7  V,  decrease  slowly 
along  the  channel.  With  the  decrease  of  cell  voltage,  current 
density  increases  quite  much  in  the  inlet  region,  and  decreases 


- 1 - 1 - 1 - 1 - ■ - 1 - ■ - 1 - 1 - 

0.0  0.2  0.4  0.6  0.8  1.0 

Specific  channel  length  (x/L) 


Fig.  3.  Current  density  and  cathode  unflooded  pore  fraction  (subplot)  along  the 
flow  channel  at  different  discharge  voltage  in  base  case.  Operation  conditions 
are  listed  in  Table  2. 


Fig.  4.  Threshold  pore  diameter  along  the  flow  channel  at  different  discharge 
voltage  in  base  case. 


along  the  channel.  Greater  current  density  gradient  is  shown 
for  lower  voltage.  This  current  density  distribution  behavior  is 
similar  to  the  experimental  results  in  our  previous  report  [29] 
and  other  literature  [30].  In  Fig.  3,  the  subplot  is  the  cathode 
unflooded  pore  fraction  at  different  cell  voltage.  It  shows  that 
the  cathode  unflooded  pore  fraction  is  smaller  than  0.7,  means 
that  all  hydrophilic  pores  are  flooded  (hydrophilic  pores  takes 
30%  of  the  pore  volume  as  listed  in  Table  2).  Along  the  chan¬ 
nel,  with  the  increase  of  water  activity  (see  Fig.  5),  smaller 
and  smaller  hydrophobic  pores  will  be  flooded  (see  Fig.  4), 
results  in  lower  and  lower  unflooded  cathode  GDL  pore  frac¬ 
tion  and  active  reaction  area,  so  current  density  decrease  along 
the  channel.  At  lower  cell  voltage,  more  water  is  produced  and 
the  unflooded  pore  fraction  begins  to  drop  at  a  nearer  place  to  the 
inlet. 

Fig.  4  shows  the  threshold  pore  diameter  for  both  anode  and 
cathode.  With  the  normal  distribution  shown  in  Fig.  2,  it  is  clear 
that  the  pore  volume  fraction  is  almost  zero  when  pore  diameter 
is  greater  than  16  |xm,  however,  for  the  calculation,  a  wide  range 
of  pore  diameter  is  taken  in  consideration  even  though  the  pore 
volume  is  very  small.  The  picks  in  the  curves  in  Fig.  4  reveal 
at  where  the  gas  is  just  in  saturation,  and  where  the  hydrophilic 
pores  are  totally  flooded.  Before  the  picks,  the  threshold  diam¬ 
eters  are  for  hydrophilic  pores,  and  the  values  after  the  picks 
are  for  hydrophobic.  In  this  figure,  the  cathode  threshold  diam¬ 
eters  fall  into  the  macropore  range,  which  influence  the  current 
density  greatly. 

Fig.  5  shows  the  species  distribution  inside  the  GDLs  and 
membrane.  It  shows  that  even  the  feed  air  is  not  saturated,  in 
the  GDL,  the  gas  is  oversaturated  from  the  inlet;  results  in  water 
activity  in  cathode  GDL  higher  than  1.0  and  membrane  water 
content  higher  than  14  on  the  cathode  side.  Oxygen  concentra¬ 
tion  gradient  in  the  inlet  region  is  higher  than  the  outlet  region, 
corresponding  to  higher  current  density  in  the  inlet  region.  On 
the  anode  side,  with  the  consumption  of  hydrogen,  water  activity 
increases  gradually  from  the  inlet  to  the  outlet.  Near  the  outlet, 
hydrogen  is  in  oversaturation. 
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Oxygen  molar  fraction  in  cathode  GDL 
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Fig.  5.  Species  distribution  in  GDLs  and  membrane  at  0.5  V  discharge  in  base 
case. 


4.2.  Cell  temperature 

Operation  temperature  impacts  fuel  cell  performance  obvi¬ 
ously.  A  higher  temperature  will  not  only  activate  the  reactions 
on  the  catalyst  surface,  but  also  accelerate  the  species  transporta¬ 
tion.  Performances  of  the  fuel  cell  at  four  temperatures  40, 50, 60 
and  70  °C  are  measured  experimentally  and  calculated  with  the 
model.  Polarization  curves  are  shown  in  Fig.  6,  where  plot  (a)  is 
the  experimental  results  and  (b)  is  the  model  results.  Reactants 
are  humidified  at  the  same  temperature  to  the  fuel  cell  reaction. 
In  the  model,  relative  humidity  of  the  reactants  is  fixed  to  be 
90%.  It  illustrated  that  the  experimental  and  calculated  polar¬ 
ization  curves  are  very  similar.  A  better  performance  is  shown 
at  higher  temperature:  higher  cell  voltage  at  a  specified  current 
density  value,  higher  limiting  current  density,  and  lower  slope 
of  the  curve  in  moderate  voltage  range  which  indicates  lower 
resistance,  which  could  be  seen  from  membrane  ohmic  resis¬ 
tivity  and  overpotential  shown  in  Fig.  7.  Relationship  between 
membrane  resistivity  and  cell  temperature  is  clearly  shown  in 
Eq.  (25).  In  Fig.  6(b),  the  cathode  activation  overpotential  is 
also  shown.  The  cathode  overpotential  difference  between  two 
temperatures  is  shown  to  be  higher  than  the  cell  voltage  differ¬ 
ence,  this  is  because  at  higher  temperature,  thermodynamic  V()c 
is  lower  from  Eq.  (28).  Voc  will  decrease  about  10  mV  when 
temperature  increases  10  K. 

4.3.  Reactants  humidity 

Because  the  conductivity  of  perfluorinated  proton  conduct¬ 
ing  membrane  relies  greatly  on  the  water  content,  reactants 
need  to  be  humidified  to  enhance  performance  of  fuel  cell.  In 
Fig.  8,  experimental  results  and  model  results  of  the  fuel  cell 
with  humidification  temperature  50,  60  and  70  °C  are  shown. 
We  assume  that  the  relative  humidity  of  reactants  out  of  the 
humidifier  is  90%  to  the  humidification  temperature.  Fuel  cell 
temperature  is  fixed  to  be  70  °C.  Use  Eq.  (12),  we  can  calcu¬ 
late  that  the  reactants  with  60 °C  humidification  is  76.8%  RH 


Fig.  6.  Polarization  curves  of  the  fuel  cell  at  different  temperature  (a)  experi¬ 
mental  results;  (b)  model  results. 


Fig.  7.  Membrane  resistivity  and  ohmic  overpotentials  at  different  operation 
temperatures. 


Membrane  resistivity  (Q*cm  ) 
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Fig.  8.  Polarization  curves  of  the  fuel  cell  with  different  reactants  relative  humid¬ 
ity  (a)  experimental  results;  (b)  model  results. 

and  50  °C  humidification  to  be  57.6%  RH  for  70  °C  reaction.  In 
Fig.  8(a),  the  experimental  polarization  curves  show  that,  with 
lower  humidification  temperature,  the  fuel  cell  behaves  worse 
performance  in  high  voltage  range  0.5-1  V,  but  has  higher  lim¬ 
iting  current  density.  Similar  behavior  is  shown  in  the  model 
results  in  Fig.  8(b)  give.  Nevertheless,  the  model  is  not  sensi¬ 
tive  enough  to  reveal  experimental  details.  For  the  experimental 
curves,  there  are  great  differences  in  high  voltage  range  (>0.8  V), 
which  is  not  shown  in  the  model  curves.  This  indicates  that 
the  catalyst  activity  is  greatly  reduced  with  lower  temperature 
humidification,  which  might  be  attributed  to  active  catalyst  sur¬ 
face  loss:  with  low  humidity,  the  proton  conductive  ionomer 
agglomerates  in  the  catalyst  layer  will  dehydrate  and  at  least 
partially  loss  their  conductivity,  which  will  make  the  platinum 
catalyst  particles  attached  on  them  insufficiently  utilized;  as  a 
result,  effective  catalyst  surface  is  reduced  and  greater  activation 
overpotential  is  shown.  Unfortunately,  this  effect  is  not  taken  into 
account  in  the  present  model  as  the  catalyst  layer  is  simplified 
to  be  an  ultra-thin  layer. 

Fig.  9  shows  current  density  and  membrane  resistivity  distri¬ 
bution  along  the  flow  channel  at  0.5  V  discharge  with  different 
reactants  humidity.  With  lower  humidity,  the  membrane  resis¬ 


Fig.  9.  Current  density  and  membrane  resistivity  distribution  along  the  flow 
channel  at  different  humidification  temperature  at  0.5  V  discharge. 


tivity  is  much  greater  at  the  inlet  region.  With  water  produced 
along  the  channel,  membrane  resistivity  decreases  and  the  cur¬ 
rent  density  increases.  When  the  gases  are  in  oversaturation, 
parital  flooding  will  make  the  current  density  drop  again.  As  a 
result,  the  current  density  distribution  for  64.7%  RH  and  76.8% 
RH  forms  a  “hill”  like  shape.  Fig.  10  shows  the  current  density 
distribution  with  64.7%  RH  reactions  feed  at  different  discharge 
voltage.  It  illustrated  that,  at  0.7  V  discharge,  different  to  the  cur¬ 
rent  density  distribution  at  base  case,  current  density  increases 
along  the  flow  channel.  With  the  decrease  of  cell  voltage,  current 
density  increase  rate  in  the  middle  of  the  channel  is  higher  than 
the  inlet  region,  and  current  density  in  the  outlet  region  begins 
to  drop  because  of  flooding.  As  a  result,  the  membrane  with¬ 
out  fully  saturation  and  GDL  flooding  effect  are  shown.  This 
current  density  distribution  trend  is  similar  to  the  experimental 
results  in  our  former  report  [29].  Although  the  operation  con¬ 
ditions  here  are  different  to  those  in  [29],  and  it  is  not  logical 
to  say  the  model  current  density  distribution  is  “validated”  by 
experimental,  the  distribution  trend  is  similar  to  some  extent.  It 
seems  that  the  dehydration  effect  is  underestimated  and  flooding 


Fig.  1 0.  Current  density  distribution  along  the  flow  channel  at  different  discharge 
voltage  with  64.7%  RH  reactants  feed. 
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(b)  Current  density  (A/cm2) 


Fig.  1 1 .  Polarization  curves  of  the  fuel  cell  with  different  air  flow  rates  (a) 
experimental  results;  (b)  model  results. 

effect  is  overestimated  compared  with  the  experimental  results. 
We  believe  that  the  this  model  will  be  more  delicate  to  predict 
the  current  density  distribution  if  it  more  details  inside  GDL  and 
catalyst  were  taken  into  account. 

4.4.  Reactants  flow  rate 
4.4.1.  Vary  airflow  rate 

If  a  fuel  cell  is  operating  with  higher  air  flow  rate,  air  stoi¬ 
chiometry  in  another  word,  not  only  oxygen  gradient  along  flow 
channel  and  through  GDL  will  be  reduced,  but  also  flooding 
will  be  lightened,  as  a  result,  a  better  fuel  cell  performance  is 
achieved.  However,  fuel  cell  efficiency  will  decrease  because 
more  energy  is  needed  for  the  air  compressor  or  air  fan.  The 
influence  of  air  flow  rate  to  fuel  cell  performance  is  shown  in 
Fig.  11(a)  and  (b)  in  experimental  and  in  model,  respectively. 
A  higher  limiting  current  density  is  exhibited  with  greater  air 
flow  rate.  Experimental  results  show  limiting  current  density 
increases  from  0.9  to  1.19  A  cm-2  when  air  flux  increases  from 
0.2  to  0.5  SL  min- 1 .  However,  the  model  results  show  greater 


Specific  channel  length 


Fig.  12.  Species  distribution  in  GDLs  and  membrane  at  0.5  V  with  0.4  SL  min  1 
air  flow  rate. 

current  density  increase.  One  of  the  reasons  for  this  might  be 
the  ignorance  of  the  current  density  difference  in  the  cross  chan¬ 
nel  direction  in  2D  model.  As  shown  in  [28],  current  density 
shows  much  difference  between  the  place  beneath  the  channel 
and  beneath  the  rib,  especially  at  high  overpotential  operation. 
The  omission  of  the  collector  will  overestimate  mass  transfer 
and  chemical  reaction  rate  beneath  the  rib  at  high  current  den¬ 
sity  operation.  The  difference  between  2D  model  and  3D  model 
was  also  clearly  shown  in  [31],  where  much  current  density  dif¬ 
ference  is  shown  at  low  cell  voltage. 

As  shown  in  the  subplot  of  Fig.  1 1 ,  where  current  density  for 
different  air  flow  rate  at  0.5  V  discharge  is  shown,  with  higher 
air  flow  rate,  cathode  flooding  is  lighter  and  current  density  in 
the  outlet  region  is  higher.  Fig.  12  shows  the  species  distribution 
in  GDLs  and  membrane  for  0.4  SL  min-1  air  flow  rate  at  0.5  V 
discharge.  Compared  with  the  base  case  in  Fig.  5,  it  is  clear 
that  with  higher  air  flow  rate,  higher  oxygen  and  water  gradient 
is  shown  in  GDL  and  membrane  in  y  direction,  and  smaller 
water  activity  in  cathode  GDL  will  result  in  lighter  flooding. 
Thus,  higher  air  flow  rate  will  be  helpful  for  water  removal  and 
flooding  prevention. 

4.4.2.  Vary  hydrogen  flow  rate 

For  the  case  of  different  hydrogen  flow  rate,  the  experi¬ 
mental  results  (Fig.  13(a))  show  almost  overlapped  polarization 
curves.  Different  to  oxygen  in  cathode  gas,  in  which  the  reacting 
agent,  oxygen,  is  less  than  20%  because  the  dilution  of  nitrogen, 
hydrogen  in  the  anode  gas  is  higher  than  80%.  Stoichiometry 
does  not  have  great  impact  on  the  concentration  of  hydrogen. 
More  importantly,  hydrogen  oxidation  reaction  is  very  rapid 
compared  with  oxygen  reduction.  Furthermore,  hydrogen  dif- 
fusivity  is  higher  than  oxygen.  These  are  the  reasons  for  very 
little  impact  shown  for  different  hydrogen  flow  rate.  The  model 
results  (Fig.  13(b))  also  show  very  small  difference  in  the  polar¬ 
ization  curves.  The  species  distribution  plot  (not  shown)  with 
0.2  SL  min-1  H2  flow  rate  at  0.5  V  shows  very  little  difference 
to  the  base  case  in  Fig.  5. 
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Fig.  13.  Polarization  curves  of  the  fuel  cell  with  different  hydrogen  flow  rates 
(a)  experimental  results;  (b)  model  results. 

4.5.  GDL  property 

4.5.1.  Average  pore  diameter 

Fig.  14  shows  the  influence  of  GDL  average  pore  diameter 
on  fuel  cell  performance,  and  the  subplot  is  current  density  at 
0.5  V  discharge.  In  Eq.  (1),  it  is  clear  that  a  smaller  average  pore 
diameter  will  result  in  higher  overpressure  for  water  condensa¬ 
tion  in  the  hydrophobic  pores.  Fig.  14  shows  performance  of  five 
average  pore  diameter  values:  6,  8,  10,  12  and  14  |xm.  Differ¬ 
ence  is  shown  in  low  voltage  range:  for  a  smaller  GDL  average 
pore  diameter,  overpressure  for  the  hydrophobic  pores  will  be 
higher,  and  that  the  flooding  is  lighter,  so  higher  current  density 
is  shown. 

4.5.2.  Pore  size  distribution 

As  normal  distribution  function  is  utilized  in  the  model  to 
simulate  GDL  pore  distribution,  the  coefficient  at  in  the  func¬ 
tion  is  modified  to  study  the  influence  of  pore  size  distribution 
to  cell  performance.  Fig.  15  shows  the  pore  size  distribution 
with  different  at  values  and  the  accumulated  volume  fractions 
integrated  from  macro  pores  to  micro  pores.  With  a  greater  at 


0.0  0.4  0.8  1.2  1.6 

(b)  Current  density  (A/cm2) 

Fig.  15.  Influence  of  GDL  pore  distribution  on  the  fuel  cell  performance  (a) 
pore  size  distribution  with  different  normal  distribution  function  coefficients; 
(b)  polarization  curves. 
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value,  the  GDL  pores  distribute  in  a  broader  diameter  range. 
The  polarization  curves  in  Fig.  15(b)  show  difference  in  con¬ 
centration  polarization  range:  with  a  greater  a>  value,  the  cell 
voltage  drops  flatter.  The  trend  is  also  shown  in  the  current  den¬ 
sity  distribution  curves  in  0.5  V  discharge  in  the  subplot.  Current 
density  begins  to  drop  at  a  nearer  place  from  the  inlet  end  for 
a  greater  a>  value,  but  higher  to  the  outlet  end.  A  broader  pore 
size  distribution  is  favorable  for  more  uniform  current  density 
distribution. 

It  must  be  acknowledged  that  the  “validation”  of  the  model 
was  limited  in  comparing  a  set  of  fuel  cell  polarization  curves 
with  experimental  results,  which  might  not  be  rigorous  enough, 
because  the  actual  physical  processes  are  much  more  complicate 
than  the  present  understanding  and  the  polarization  curve  is  only 
a  gross  behavior  of  all  the  microcosmic  processes.  Although 
the  current  density  distribution  trend  is  similar  to  experimental 
results  to  some  extend,  the  detailed  behavior  predicted  by  the 
model  has  not  been  validated  by  experiments,  and  as  such,  the 
predictions  do  not  necessarily  represent  the  actual  behavior  of 
an  operating  fuel  cell. 

5.  Conclusions 

A  two  dimensional  partial  flooding  model,  in  which  GDL 
pore  size  distribution  is  first  taken  into  consideration,  is  devel¬ 
oped  to  study  the  influence  of  liquid  water  to  the  performance 
of  PEMFC.  In  this  model,  liquid  water  produced  is  considered 
to  condense  first  in  hydrophilic  pores  and  then  in  hydrophobic 
pores  if  water  vapor  pressure  is  higher  than  the  condensation 
pressure  for  the  pores.  The  partial  flooding  will  reduce  the 
active  reaction  area  in  turn.  Use  this  model,  different  operation 
conditions  are  studied  and  the  model  results  show  reasonable 
accordance  with  experimental  results. 
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